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We propose an efficient method to compute reaction rate constants of thermally activated pro- 
cesses occurring in many-body systems at finite temperature. The method consists in two steps: 
first, trajectories are sampled using a transition path sampling (TPS) algorithm supplemented with 
a local Lyapunov bias favoring diverging trajectories. This enhances the probability of observing 
rare reactive trajectories between stable states during relatively short simulations. Secondly, reac- 
tion constants are eventually estimated from the unbiased fraction of reactive trajectories, yielded 
by an appropriate statistical data analysis tool, the multistate Bennett acceptance ratio (MBAR) 
package. 

In order to test our algorithm, we compute reaction constants for structural transitions in LJ38, 
a well studied Lennard- Jones cluster, comparing our results to values previously reported in the 
literature. Additionally, we apply our method to the calculation of reaction rates for the migration 
of vacancies in an a-Iron crystal, for temperatures ranging from 300 K to 850 K. Vacancy diffusion 
rates associated with activation barriers at finite temperature are then evaluated, and shown to be 
substantially different from values obtained using the standard harmonic approximation. 

I. INTRODUCTION 

Rare events are physical events occurring with low probability in numerous many-body systems, encompassing 
a wide range of phenomena, from condensed matter physics to proteins: they often involve thermally activated 
processes, i.e. passages between different stable configurations of the system that are separated by free energy 
barriers. Evaluating the frequency of these rare events from numerical simulations can be very time consuming, as 
the probability of observing one of these passages is very low: for instance, the migration of a vacancy in a-Iron at 
500K typically happens every microsecond, while the usual time steps for molecular dynamics is of the order of a few 
femtoseconds. In the last decades, several techniques have been developed in order to accelerate the dynamics, and 
enhance the probability of observing infrequent events during a short simulation, based on importance sampling. ^iliIS 

Among these techniques, transition path sampling 3 (TPS) allows to estimate the frequency of rare events by means 
of path ensemble averages of physical observables: reaction rates for appropriate time scales are indeed evaluated from 
the ratio of the number of reactive trajectories on the total amount of paths sampled. However, a sufficient number 
of reactive paths has to be observed, in order to obtain a reliable statistics. 

Herein, we propose a transition path sampling algorithm where the fraction of reactive paths sampled is enhanced 
using an adequate bias that favors diverging trajectories. It is indeed possible to show 8,33 that reactive paths we want 
to sample share important features with diverging trajectories observed in chaotic systems. Therefore, a suitable 
parameter that quantifies the divergence of dynamic trajectories can be exploited also to bias a path sampling 
algorithm aimed to reproduce reactive paths. 

The main instrument proposed in the literature to quantify chaotic behavior of dynamical systems is the evaluation 
of Lyapunov exponents, 34 that are usually employed to estimate the sensitivity of deterministic systems to small 
changes in initial conditions. For this feature, they have been widely studied,— both analytically and numerically, 
in hamiltonian as well as in nonlinear systems of small dimensions. Moreover, the use of Lyapunov exponents to 
characterize numerically phase transitions in finite size systems has been extensively explored in the past years, and 
many noticeable results have been obtained in the early 90'si^ ~ 41 ' 60 

Resorting to Lyapunov exponents in order to achieve numerical ergodicity and localize saddles and transition paths 
has been done recently with the Lyapunov-weighted dynamics method proposed by Tailleur and Kurchans^ in this 
sampling scheme, a set of clones are copied or deleted depending on a probability weight computed from quantities 
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related to Lyapunov exponents. After this work, the paper of Geiger and Dellago^ 3 . have shown how to couple the 
chaoticity features of a dynamical system to a TPS technique for sampling deterministic trajectories, using an indicator 
for diverging trajectories borrowed from studies on planetary systems^ 3 -, the relative Lyapunov indicator (RLI). 

In this paper, we present a chaoticity indicator different from RLI, and based on local Lyapunov numbers, that 
are quantities closely related to Lyapunov exponents. This indicator is used to introduce a bias in the path sampling 
scheme, thus obtaining a Lyapunov biased TPS (LyTPS) method that will be in the sequel applied to complex 
many-body systems, like the well known optimization benchmark model LJ38. Furthermore, we show how reaction 
rate constants can be recovered from biased TPS quantities, resorting to an appropriate statistical analysis to unbias 
reaction rate values computed in a LyTPS framework. 

This paper is organized as follows. In Sec.Qll we first recall the basic concepts of Lyapunov exponents for dynamical 
systems. We will then briefly review the use that has been made of them in numerical algorithms to characterize 
phase transitions, or in importance sampling contexts. We then expose how to use local Lyapunov numbers in the 
context of a Transition Path Sampling to determine saddle points and reactive paths (Sec. IIIip . Reaction constants 
are computed from the fraction of reactive paths using the Bennett-Chandler approach 52 of population correlation 
functions and the standard TPS technique* 2 -; we also explain how unbiased reaction constants are recovered from a 
Lyapunov biased algorithm thanks to the multistate Bennett acceptance ratio^ (MBAR) method (Sec. IIVj) . Finally, 
numerical results concerning the application of our method to solid-solid structural transition in LJ38 and vacancy 
migration in a- Iron are presented (Sec. |V|) . 

II. LYAPUNOV EXPONENTS IN DYNAMICAL SYSTEMS 

We briefly present here the theory of Lyapunov exponents, mainly following Ott 3 — ; then we propose a formulation 
allowing the use of these exponents in importance sampling techniques. 

A. Discrete dynamics and numerical applications: state of the art 

In numerical applications, dynamics are discrete: the evolution of state vector x at time step n is described by a 
mapping x„ + i = M(x„), where M is a matrix expressing the system evolution from one time step to the following. 

Let the system be at time t = in an initial position xo , and let (5xo be a small perturbation applied to this initial 
state. The dynamics of such a perturbed system can then be denoted using a new state vector x = x + <5x, whose 
time evolution will be 

x„ + i = x n+ i + 5x„ + i = M(x n ) = M(x„ + (5x„) (1) 

For a sufficiently small perturbation <5x, it is possible to linearize the map M as 

M(x„ + <5x n ) = M(x„) + Z?M(x„) • <5x n + 0(<5x 2 ) (2) 

The time evolution of a small perturbation to the initial state vector then reads 34 

6x n+ i = DM(x„) • 6x n (3) 

where Z?M(x„) is the Jacobian matrix of the map M. Inserting particular solutions <5x„ = e[A]" in Eq. ([3]), one finds 
an eigenvalue equation 

DM(x n ) ■ e = A ■ e. (4) 



In the discrete case, the eigenvalues A^ of Z?M(x„), solutions of Eq. are called Lyapunov numbers rather the 
Lyapunov exponents, and trajectories are unstable for |Afc| > I, and stable otherwise. The largest eigenvalue A MAX 
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will be associated to the eigenvector e MAX indicating the direction of maximal growth of the perturbation i5x. We 
introduce the matricial product 

M» = M(x„_ 1 )-MW (5) 

between the Jacobian matrices of the hamiltonian map at successive time steps, and we express the perturbation at 
time step n with respect to the initial perturbation <5xo as 

<5x„ = J DM"(x )- ( 5xo (6) 

Defining with ||<5xo|| the Euclidean norm of <5xo in phase space, the Lyapunov exponents h, given the initial condition 
Xo and the initial perturbation orientation Uo = <5xo/ ||<5xo|| are 

h(xo, u ) = lim - In (^4) = Iim - In \\DM n (x ) • u || (7) 

n~¥oo n \\\OXq\\ J n-Joon 

For A-body hamiltonian systems, having 3N degrees of freedom, x(t) is a ^-dimensional state vector accounting 
for both the positions and momenta of the N particles. For the ^-dimensional hamiltonian map there will be 6N 
Lyapunov exponents, usually ordered in literature from the largest to the smallest (hi > ■ ■ ■ > /i6iv)- As stated in 
Ref. 34 , Oseledec's multiplicative ergodic theorem 4 " guarantees the existence of the limits used in the definition of the 
Lyapunov exponents in Eq. under very general circumstances. 

The Lyapunov exponents are related to the aforementioned Lyapunov numbers as = exp [hk]- From Eq. ([?])■ we 
can also define finite-time Lyapunov exponents: 

Mx , uo) = - In (^4) = - In ||£>M»(x ) • u || (8) 
n V||dx ||/ n 

For long enough times, the greatest Lyapunov number Ai will give the dominant contribution to the perturbation 
evolution, and the associated eigenvector ei indicates the direction of maximum growth of the perturbation <5x. 

We stress here that finite-time Lyapunov exponents are calculated for a given Xo and that, strictly speaking, 
their values do depend on the initial orientation Uo. It is however shown that the largest exponent fti(xo,Uo) is 
approximately independent of the choice of Uo in Hamiltonian ergodic system o 27 ' 34 , while the complete spectra of 
finite-time exponents can be determined using specific numerical techniques. 47 

In numerical simulations, only finite-time Lyapunov exponents can be estimated, due to limited amount of CPU 
time. To evaluate h n from Eq. ([5J we should compute the matricial product of Eq. ([3]). For systems with many 
degrees of freedom, a calculation of this matrix product is not possible analytically, and is numerically expensive. The 
solution most followed in literature consists in directly evaluating the quantity 

ft R (x 0) u ) = -lnf{^4 > ) (9) 



n \P x o|| 

given by the distance ||5x n || between two nearby dynamical trajectories, the first one started from the initial state Xo 
and the second one from the perturbed configuration xo = xo + <5xo after n time steps. 

The use of Eq. ([§]) as a mean to evaluate finite-time Lyapunov exponents has two main drawbacks: the need of 
computing two trajectories to evaluate a single Lyapunov exponent, thus doubling computational cost, and the fact 
that values obtained for h can be sensitive to initial conditions, because of the dependence of the computed finite-time 
Lyapunov exponents from the choice of the orientation of initial perturbation Uq, as recalled above. 

Several numerical strategies have been proposed to bypass these two problems. The tangent space method^ ' 16 ! 47 
assigns to each state x t of the trajectory started in Xo a vector u(x t ). These vectors are computed from the local 
hessian matrix of the hamiltonian mapping, and their norms indicate the distance between the current trajectory and 
the perturbed one, i.e. u(i) ~ 6x(t). As these distances evolve exponentially (see Eq. (|7J|), the lengths of the vectors 
u can quickly diverge or vanish: a reorthonormalization of u at each time step is therefore required, for instance 
with a Gram-Smith algorithm. This method has been implemented in the literature^, for instance in the context of 
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Lyapunov weighted dynamics** - — To make this algorithm independent of the choice of the first vector u(xo), one could 
integrate the equations of motion backward in time from x for a duration r, and then reintegrate the evolution of 
u(x t ) forward until t = 0(see Ref£). In this way, u(x ) would be automatically oriented in the direction of maximum 
growth. However, the duration r should be long enough to ensure the loss of correlation between the orientation of 
u(x_ T ) and u(xo), thus requiring the computation of long trajectories at sustained computational cost^^ 3 - 

Another evaluation of finite-time Lyapunov exponents can be obtained with the Relative Lyapunov Indicator (RLI) , 
elaborated by Sandor et o/.— in the context of planetary trajectories, and further used in a Lyapunov weighted path 
sampling scheme^ 3 .. The main idea is to compare finite-time Lyapunov exponents h for trajectories starting very close, 
say in xo and Xo + Axo . The difference between finite-time exponents at time step n can be written as 

Aft n (x ,u ) = |M x o + Ax ,u ) - h n (x ,u )\ (10) 

and will in general undergo strong fluctuation, 43 which can be smoothed by an average over AT trajectory steps: in 
this way one defines the RLI as the quantity 

1 M 

R(Xq, U ) = jj. ^2 A M x o, u ) (11) 
i— 1 

This average over the entire trajectory length reduces^ 3 - the dependence of the computed finite-time Lyapunov expo- 
nents on the orientation of initial perturbation, but introduces an additional dependence on Ax . Both finite-time 
Lyapunov exponents are calculated evaluating the distance between two trajectories evolving close to each other (in- 
stead of the matricial product of Eq. ([5])): in terms of computational cost, four trajectories are computed to obtain a 
single RLI. 

In the following, we propose a faster and orientation-independent way to evaluate the divergence of hamiltonian 
trajectories, alternative to tangent space method and RLI, to be used in the path-sampling scheme described in 
Sec. HIl 



B. Hamiltonian dynamics 

We restrict our focus to systems with deterministic dynamics governed by an hamiltonian T~L = 2m ^ ^"(l) - 

The time evolution of the state vector x = (q, p) (also indicated as "hamiltonian flow"—) directly follows from 
Hamilton equations, 

/ r> t \ / d?*(q»p) \ 

(12) 




dp 

The evolution of a small perturbation of positions and momenta for a standard hamiltonian is then given by 



^ x = I a 2 v(q) Q ] 5x ( 13 ) 

where G is the 3iV x 37V inverse mass matrix of elements Gij — Sij/rrii. 

Let us discretize the hamiltonian dynamics in Eq. (|12[) with the velocity Verlet algorithm: 

v , 7) 1 dt 

dt 

%n+l = Qi,n H Pt,n+l/2 (14) 

m 



1 dV(q n+1 ) 

Pi,n+1 — Pi,n+l/2 ~ 



2 dq in 



+ 1 
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This algorithm is accurate to second order and numerically stable^. It will be used in Sec. IHII to generate dynamical 
trajectories in numerical applications. The perturbation (5x n+1 can be evaluated with respect to <5x„ using Eq. ([3]), 
where the jacobian matrix of the hamiltonian mapping for a velocity Verlet discretization reads*^ 



£M(x„) 



I - £rH(x n ) 



dt 
2m 



{H(x„) + H(x n+1 ) [i - ^H(x n )] } I 



dt A 
2m 



G 

H(x„ +1 ) 



(15) 



In the upper right and bottom left blocks we introduced the hessian matrix of the potential energy H at states x ra 
and x„ + i, respectively. 

The jacobian matrix of Eq. (fTS")) obtained with the velocity Verlet scheme of Eq. ([TJJ contains the hessian matrices 
at steps n and n + 1. Therefore, manipulating DM(x„) is numerically expensive. A simpler expression for the jacobian 
matrix can be obtained from the less accurate Euler discretization algorithm: Eq. (TT2"j) becomes a set of 6N coupled 
equations of motion 



Qi,n+1 — Qi.n + dt ■ 
Pi, n+1 = Pi,n - dt ■ 



mi 
dV 
dqi, n 



where i = 1, 3N, so that the jacobian matrix of the hamiltonian map reads 



>-6N 



dt 







G 



-H(x„) 



(16) 



(17) 



Again, the perturbation 5x„ at each time step can be evaluated by inserting Eq. p7|) in Eq. ([3]). 

The difference between the jacobian matrix Z?M(x„) of Eq. (|17l) and the one of Eq. (fTS")) consists in second order 
terms. However, it is numerically less expensive to manipulate the former than the latter, as _DM(x„) of Eq. (fTT|) 
requires to evaluate the hessian only at time step n. In the following, we will be interested in computing the eigenvalues 
of DM(x„), in order to obtain a bias favoring reactive trajectories: this bias will be removed at the end, so it would be 
useless to spend CPU time to accurately evaluate the jacobian matrix. Therefore, accordingly to Ref. 4 — , we consider 
the Euler scheme (Eq. (|16|) ) precise enough for our purposes, and we use the first order Euler discretization of Eq. \YJ\ 
to compute DM(x n ). 



C. Maximum local Lyapunov numbers 

Using the discretized hamiltonian dynamics given in Eq. (|16j) . we proceed by computing at each time step the 
maximum local Lyapunov number—, given by the largest eigenvalue of DM(x„) (Eq. (|17I) V The 6N eigenvalues A„ 
of BM(x„), computed at time step n, can be obtained writing Eq. Q as a secular equation 

P(A„) = det {A n I 6N - M(x„)} = 0. (18) 

whose solutions are 57V pairs of eigenvalues A„ of Z?M(x„) , because of the simplectic properties of the hamiltonian 
mapping matrix Mi- These eigenvalues are given by the expression 40 



Af, n = 1 ± dty/-X j>n Vj = l,...,3N (19) 
where the \j, n correspond to the eigenvalues of the mass- weighted Hessian H'(x„) of components 

tt' / \ 1 d 2 V(x n ) 

H ij(x n J = - — . (20) 

v /m7ro~ dqidqj 

Eq. (fT9"|) shows that at each time step n the jacobian eigenvalues A 3 - n , i.e. the local Lyapunov numbers, depend on 
the potential energy surface through the hessian eigenvalues Aj. n : unstable configurations x„, such as saddle points are 
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characterized by the occurrence of some negative Xj <n , and correspond to real and positive local Lyapunov numbers 
Aj n . Conversely, stable states have all Aj.„ positive or imaginary local Lyapunov numbers with unitary real part. In 
the following path sampling scheme of Sec. Mil we neglect the imaginary part of Aj n given by stable states. This is not 
an issue because it is sufficient to determine real and positive Lyapunov exponents to characterize unstable dynamics 
which we are interested in, and the imaginary part of the jacobian eigenvalues Aj iU can be neglected. Indeed, the 
imaginary part of the Lyapunov exponents does not affect the stability of the system, but only indicates if dx n is 
spiraling clockwise or counterclockwise. 

At each time step n, the most negative eigenvalue of the hessian matrix A™ 71 gives the eigenvalue of the Jacobian 
matrix _DM(x„) with the largest real part. We write the maximum local Lyapunov number as 

a max = 1 + dtv / max ( 0; -A™ n ) (21) 

such that Eq. (|2Tj) entails A^ = 1 for stable configurations, where all A are positive, and A^ AX > 1 for unstable 
configurations having a negative spectra. 

Eigenvalues A„ are found by diagonalization of the hessian matrix H. This can be computationally very expensive 
for systems with a large number of degrees of freedom. One efficient solution consists in extracting only the lowest 
eigenvalue A™ m needed to evaluate A^ AX using the Lanczos algorithm^. This iterative algorithm finds extremal 
eigenvalues of any matrix with a reduced computational cost, diagonalizing only a submatrix of the initial one (see 
for example appendix A of Ref.— for details). As pointed out in Ref. — , a 15 x 15 Lanczos submatrix is sufficient to 
detect negative eigenvalues. Moreover, it is possible to decrease the submatrix size to as little as 4 x 4 by verifying 
at each iteration that the Lanczos solution is stable; if not, repeat the calculation until a the solution is converged. 
Hence, the most negative eigenvalue, corresponding to the most unstable direction of the potential energy surface at 
a give system position in the phase-space at a given instant can be extracted in a very few iterations. This is the 
numerical method we will apply in the following to evaluate A^f AX . 



D. Lyapunov indicator for dynamical trajectories 

A dynamical trajectory is defined as an ordered sequence of states in phase space separated by a small time 
increment St, and denoted as z = {xq, x T }, i.e. a path of total length r composed by Af = £ + 1 state vectors. We 
introduce a Lyapunov indicator for path z 

L(z) = ±lnf[A™ AX (22) 

n— 1 

given by the average of the maximum Lyapunov number of Eq (|21[) over the whole trajectory. 

The verification of the difference between finite-time Lyapunov exponents estimated by RLI or tangent space method 
and L(z) from Eq. (|2"2")l is beyond the scope of this article. We stress instead that the indicator proposed in Eq. (|22l) . 
being based on the hessian spectra A, is strictly related to the topological properties of the potential energy surface, 
thus encodes local information on the stable or unstable configurations sampled in phase space by a given trajectory. 
Henceforth, we consider this Lyapunov indicator suitable for importance sampling techniques. 

The largest local Lyapunov number of Eq. (|2Tj) were evaluated by Hinde et aL— to study the dependence of the 
Kolmogorov entropy on the potential energy surface of small Lennard- Jones clusters. In that work, the Lyapunov 
exponents derived from the eigenvalues of the jacobian matrix of the hamiltonian mapping are summed over trajectories 
of different lengths, thus obtaining - using Pesin's theorem— - an estimation of Kolmogorov entropy. Using Eq. ([21]) to 
evaluate finite-time Lyapunov exponents was proven to be a quite successfull approach, as it yields enough information 
to quantify the degree of instability of phase space trajectories. This supports the idea that evaluating the global 
Lyapunov exponent from local Lyapunov numbers allow to correctly characterize the chaotic properties of the system. 
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III. TRANSITION PATH SAMPLING WITH A LYAPUNOV BIAS 



The idea of sampling the phase space of a many-body system through paths generated by molecular dynamics, using 
the Metropolis algorithm, has been developed by Dellago and coworkers . 2 i 36 i 45 This approach is called transition path 
sampling (TPS). In this present work, we introduce a bias in the TPS algorithm in order to favor the sampling of 
reactive trajectories. The bias is proportional to the Lyapunov indicator L(z) obtained in Eq. (|22[) . 

Each path z constrained to start in a reactant basin A is equipped with a probability density 

= 4 exp {aL(z) - /3H(z)} ^(x ) (23) 

where Z\ is the partition function on the biased trajectory ensemble, and function </?J|(xo) is an additional term linking 
the initial state Xo of the path to state x^. Different choices for <^ are possible, for instance <^(xo) = /ia(xo), where 
Ha is an indicator function on A, such that 

{1 x G A 
x^A 

or 

<^(x ) = exp j-^ Q (x -x A ) 2 J (25) 

that accounts for having a tunable spring of stiffness K a linking the origin of path z to state A. In this last case, the 
stiffness parameter k" can be tuned to counterbalance the strenght of the bias. We denoted in Eq. (|23|) for simplicity 

r/St-l 

exp{-/3H(z)} = p(x ) Y[ 6 [ x d+i)st - (t>st(*ist)] (26) 

i=0 

as the dynamical path probability arising from the deterministic propagation of the trajectory. In Eq. (|26[) p(jx-a) = 
exp(— /3'H(xo)) is the unnormalized canonical distribution at inverse temperature f3 from which the initial configuration 
is selected, and <fi is the temporal propagator x t = </> t (xo) associated to the deterministic dynamics 3 ^. 

In this biased ensemble, choosing positive a enhances the probability weights of trajectories with a large Lyapunov 
indicator L(z), favors via Eq. (jT§]) the sampling of reactive paths passing over saddles and unstable directions of the 
potential energy landscape. On the contrary, choosing negative a mainly restricts the sampling of non reactive or 
regular trajectories within stable basins. 

The distribution V\ is approximated by a Markov chain of M steps constructed by importance sampling, by means 
of the Metropolis algorithm. The sampling is done in the following way: at Markov chain step m, starting from the 
current path z m , a trial path z is generated with probability P gen . Then, the trial path is accepted with a probability 
Pace and added to the Markov chain as z m+1 = z; otherwise, if the trial path is rejected, z" 1+1 = z m . To ensure the 
convergence of the Markov chain towards the equilibrium distribution Va, we impose that the probability n [z — > z'] 
to transit from a path z to a different path z' satisfies the detailed balance equation 

^[z]7r[z^z']=^[z']7r[z'^z] (27) 

Taking account of the generating and acceptance probabilities P gen and P aC c, the transition probability 7r reads 

7T [Z -> z'] = P 9^ [' 2] ^ z')Pacc [ Z ~> ^ + 5 ( Z ~ Z )(! ~ P «cc [z ->■ S])} (28) 

z 

where 5 is the delta distribution and z' is either the old path z or the proposed path z. The acceptance probability 
can be constructed from Eqs. (|27[) and (|2"5f as the Metropolis acceptance 

that is widely used in numerical simulations, as it has the main advantage of maximizing P acc . 



8 



Shooting algorithm 

The standard shooting algorithm for deterministic dynamics is obtained by (i) selecting a state Xj< of the current 
trajectory, (ii) perturbing the momenta of each particle, and then generating from this selected state two segments, 
one backward of duration t' and the other one forward of duration t — t', in order to get a trial trajectory z of same 
duration, (hi) accepting or rejecting the new trajectory z. For deterministic dynamics, the total energy is constant. 
The perturbation step (ii) is done with the algorithm proposed by Stoltz 5 -, where trial momenta p are obtained from 
old momenta p as 



p = ep + v 7 ! - £ 2 6p (30) 

where e is a tunable parameter and Sp is drawn from a Gaussian distribution of variance ksT. The probability of 
having trial momenta p from old momenta p with the Stoltz proposal for shooting^ 5 - is written as 

(\ 3N ^ 
1 \ J (p-6 P ) T (p-ep) \ 
7W^) 2(1 -e») I" (31) 



and ensures the detailed balance condition 

exp{-/3H(x t >)}p(Pt> -> Pi') 



(32) 



exp{-/3"H(x t 0}p(Pt' ~> P*') 

hence the probability flux between the current and perturbed momenta is balanced and the hamiltonian distribution 
is preserved. 

The probability p ge „[x t / — > x t /] of obtaining the shooting point x t / for the trial trajectory from state Xf selected in 
the current trajectory reads 

Pgen[xt> ~* Xt'] = P (Pi' Pi') (33) 

as only momenta are perturbed at the shooting point, while positions are left unchanged. 

Inserting Eq. (|23[) in Eq. (|29p and neglecting numerical approximations in the integration from xj to xo (that indeed 
give in computations "H(xJ) 7^ "H(xo)) we have the Metropolis acceptance rule 

P acc [« h. z] = min (l, exp {aL(z) - ^(x )} ^(x ) Pge »[x t , -> x,j 1 

1 J I 'exp{aL(z)-/3^(x )}^(xo)p S en[xi' ^x t /]J 1 ' 

where terms deriving from the forward and backward generation of the current and trial path have cancelled because 
of the unit phase space compressibility of the Newtonian dynamics (see Ref<^ for details). 

Using in Eq. (|34| the property of the Stoltz proposal (Eq. ([32])) and neglecting numerical approximations in the 
integration from xj to Xoj Eq. (|2T)|) further simplifies in 

Pace [z -> z] = min{l,exp{aL(z) - aL(z)} [Va(*o)/V1( x o)]} . (35) 



Shifting algorithm 

The second Monte Carlo move in trajectory space is based on the shifting algorithm, supplemented with a waste- 
recycling estimator^ Af trial trajectories Zj are constructed from z as follows: the duration of the trajectory z is 
doubled selecting a random duration vSt and integrating two segments backward and forward, starting from Xo and 
x r , along v and M — v time steps, respectively. Adding these segments to the current trajectory, one obtains a 
"buffer" trajectory ( = {Xnst} -v< n <2fS-u 01 total duration 2J\fSt, containing J\f possible trial paths. The conditional 
probability of obtaining the "buffer" trajectory £ starting from the current trajectory z is indicated as -P C ond(Cl z ) • 
The probability weight of each trial path Zj is then 

Wi) = ^exp{aL(2 i )-j8H(2i)}^(^) ( 36 ) 
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if a constraining function ^^(xo.j) linking the initial state icoj of each trial path Zj trajectories to state A is used, as 
in Eq. ([25)1. Index j runs over the A/" possible paths on the "buffer" trajectory. 
We introduce an action 

— s a j = aL(zj) — /3H(zj) (37) 
so as to rewrite the trajectory probability in the biased ensemble (Eq. (I36|) ) as 

£a( z j) = ^^A(xoj)exp{-s QJ }. (38) 

We now define the selecting probability P se i(zj\Q of selecting a trial trajectory z,,- from the buffer trajectory £ as 

<^(x 0j )_exp{-s a j} 

n a 

where we introduced the Rosenbluth factor 

M 

K> a = ^2 ^S(xj) exp [-s QJ ] . (40) 
i=i 

Resorting to Bayes theorem, P se ; can be written as the posterior likelihood probability^ of having Zj given the 
"buffer" trajectory (: 

PselizjlO = Pc ^ff j) ^ j) (41) 
raarg \S> ) 

where P CO nd (Cl^j) 1S the conditional probability of constructing a "buffer" path £ from the trial path, and because of 
the deterministic dynamics Pcand((\zj) = l/A/\ PmargiO is the marginal probability associated with the buffer path: 
comparing Eq. (|4"T)) with Eq. ([3"9")l . we see that 

Kar g (0=^^^- (42) 

Let us now consider a Monte Carlo move between two paths both contained in the buffer trajectory £ , i.e. from 
the current path z to Zj, whose associated transition probability 7r[z — > Zj], as in Eq. (|27[) . obeys a detailed balance 
with respect to the prior distribution "P^: 

VI [Z] 7T [Z -> Z,] = ^ [%] 7T [ Zj -> Z] . (43) 

Defining the transition probability 

7T [z -> %] = P se *(%IC)Pc «d(C|z) (44) 

the detailed balance in Eq. ([43)1 can be recasted as 

P sel {z i \QP cond {C\z)V'l{z) = P sel (z v \OPcond((\Zj)pA&) (45) 

where z = z v and P; e ;(z„|£) is the probability to transit from Zj to z. As a consequence, the shifting procedure leaves 
the probability distribution V\ invariant. Moreover, recalling that the deterministic dynamics entails P CO nd(C\ z ) = 
Pcond{C\zj)i the detailed balance in Eq. ([45)) further simplifies into 

p.eifoicra*) = Psei(z„\on(^)- (46) 

We conclude this section by extending the detailed balance of Eq. (|43)) for trajectories z and Zj belonging to two 
different buffer paths ( and ( respectively. Writing as in Eq. ([4T)) expressions for P se i(z\Q and P se ;(zj|£), and using 
these results in Eq. (|4"3")l we obtain 

PZ arg (OPsei(z\On [z Zj ] = P r ? larg (C> sei (z,|C> z] (47) 

Eq. (|47l) is indeed a detailed balance condition for the alternate shooting and shifting moves, and samples the buffer 
trajectories £. This implies that the distribution P^arg i s invariant along the sampling algorithm. P^arg I s therefore 
suitable to be used as an input path probability weight required by the unbiasing algorithm MBAR, that will be 
presented and used in the following Sections. 
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IV. REACTION-RATE CONSTANTS CALCULATION 
A. Rate constants theory 

Here we briefly recall how reaction rates are computed in the TPS framework. The reactivity of the sampled paths 
is recovered by means of a time correlation function with respect to initial and final states of the paths: a trajectory 
is said to be reactive if it starts in the reactants A basin and ends in the products B basin. The time correlation 
functionii2& is 

= (WW) ( 

(Au(xo)) 

where brackets (•) indicate averages taken over the equilibrium trajectory ensemble and the indicator function hn is 
defined as in Eq. (f2"4"|). C(t) can be understood as the conditional probability of observing a trajectory of duration t 
ending in state B, knowing that it started in state A^ i 36 ' 57 , and approaches its asymptotic value exponentially as 

C(t) « p e £ (1 - exp {-t/ Trxn }) (49) 

where p e g is the equilibrium occupation probability of state B, and the parameter r rxn = {Ua—^b + &b-ka) 1 is the 
characteristic reaction time of the system, given by the forward and backward reaction constants fc^^s and ks^A, 
respectively. 

Note that the basic assumption required to compute reaction rate constants of rare events from the correlation 
function C(t) is the presence of a well separated time scale for processes occurring between 'fast' intra-funnel relaxation, 
having a typical time constant r mo ;, and activated processes indicating passages between funnels, needing a much longer 
time scale, of the order of r rxn —. 

For times in the intermediate time regime r mo i < t -C r rxn , the correlation function in Eq. (|49[) can indeed be 
expressed by its first order expansion 

C{t) « k A ^ B t (50) 

where the detailed balance condition Ua-^bPa = ^b^aPb nas been used to eliminate p e S. Hence, the slope of C(t) 
for this intermediate time regime gives direct access to reaction rates. The reactive probability flux flowing from state 
A towards B per unit time, defined by k(t) = dC Jf^ , displays a plateau corresponding to the forward phenomenological 
reaction constant kA^B^ 

Let us point out that these results, obtained with a macroscopic 'flux over population' probability approach, can be 
recovered by a Bennett-Chandler formalism—, based on microscopic quantities (positions and momenta, see Refj* 7 -). 
Finally, we recall that the phenomenological forward rate constant kA—^B computed from Eq. (|50"|) is related to the 
rate krsr given by the transition state theory (TST) 10 i 52 

k(t) = n{t)k TS T (51) 
where n(t) is the transmission factor, and fcrsr reads 

krsr = exp (-P^a^b) (52) 

with AFa->b the height of the free energy barrier separating basins A and B, h the Planck constant and /3 = 
the inverse temperature. The transmission factor k is always lower than one, and is introduced to take into account 
trajectories started in basin A that reach the saddle point but fall back to state A, instead of ending in state B: these 
occurrences are called recrossings events. The transmission coefficient usually reaches a steady value, depending on 
the temperature and the reaction coordinate chosen to localize the barrier. At this plateau value of R < 1, and the 
reactive flux corresponds to kA^B = Sfcrsr, always lower than the TST value. 
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B. Rate constants with biased sampling and waste-recycling 

In numerical experiments, the computation of reaction constants by direct evaluation of C{t) at times longer than 
the intrafunnel relaxation time r mo ; means performing very long molecular dynamics trajectories. To estimate C(t) in 
the TPS framework computing relatively short trajectories, one resorts to a factorization of the correlation function in 
(i) a static quantity related to kxsr-, thus dependent on the free-energy barrier and calculated through an umbrella- 
sampling technique, and (ii) a dynamic factor related to «(£) given by the time derivative of the probability of reaching 
basin B at times shorter than the whole trajectory length (see for instance Refs. — ^) 

We propose herein a variant strategy: reaction constants will be calculated directly by averaging indicator functions 
on short trajectories, as in Eq. (|48j) . once the fraction of reactive paths is significantly enhanced by introducing an 
appropriate bias favoring reactions between basins A and B. 

The correlation function in Eq. (l48t can be intended as an average over all performed trajectories - i.e. over all 
successive steps m of the Markov chain - of the reactivity A(z m ), defined as 

A(z m ) = Mx™)Mx;") (53) 

For an unbiased TPS simulation, the probability V\ of Eq. (|23|) with a = is sampled: we have C(t) = (A} , where 
brackets (-) correspond to averages over a canonical trajectory ensemble, and the trajectory distribution of Eq. (|23[) 
ensures (/i J 4(xo)) = 1 because of ip^. 

In a context of biased TPS, averages on Markov chains are taken on the biased trajectory distribution, hence 
we denote biased averages of the correlation function as C a (t) = (A) , where index a accounts for the current bias. 
Herein, index a will indicate all observables obtained from a biased distribution, where a = stands for the equilibrium 
canonical ensemble average. We estimate the correlation function C a {t) = (A) a using an estimator denoted by I^f , 
which consists in taking the average over the Markov chain of lenght M as 



l « M = If E A « ( 54 ) 



A(z™) = A™ denotes a reactivity value coming from a biased trajectory z a distributed according to V% Estimates 
given by l£f are however not optimal^. 

A waste-recycling (WR) estimator 70 i 71 is associated to the multiple proposal sampler for the shifting move of 
Sec. EH to obtain more accurate estimates. Waste recycling consists in including information about all possible 
proposals contained in the buffer trajectory £. The reactivity at Markov chain step m in a given ensemble a has to 
be first averaged over the A/" trajectories contained in the buffer path £™ as 



X = E 4Sj ^ X °^P[ s ^ = J2 AZ^ A (xft) exp [- s - + SI 



3 = 1 3=1 



(55) 



where we write A a — A(C™) and define the Rosenbluth factor 

-RZ = exp[-SZ] (56) 

as proportional to the marginal probability PmargiC™) of Eq. PO)) . associated to the "buffer" trajectory £™ corre- 
sponding to Markov chain step m. Correlation functions for reactive paths are therefore estimated in a way similar 
to Eq. (|5"i|) . with the WR estimator 

J M U1 _ E^ 1 Ef =1 < J ^(x^)exp[-C J +^] _lf T 
and again C a (t) ~ J a [A]. The calculation of rate constants fc^-^s f° r the given a-ensemble follows from Eq. (f50|) . 
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C. Unbiasing rate constants: the MBAR algorithm 



A suitable unbiasing algorithm is needed in order to recover canonical ensemble correlation functions from reactivity 
values witnessed in a Lyapunov biased path ensemble. The canonical equilibrium values of Co(t) can in principle be 
obtained by estimating reactivities A™, computed in any a-biased path ensemble, resorting to an adequate unbiasing 
algorithm. We define an unbiasing WR estimator Se,a, where the left subscript 9 indicates the ensembles in which we 
are interested in measuring averages, while the right subscript a refers to the ensemble that our Lyapunov biased TPS 
will effectively sample. Equilibrium values for Co(t) (9 = 0) are retrieved from reactivities computed in any a-biased 
ensemble as 

cm - C Mi - ^%^«^'<*^ bfS* + g + m 

Unbiasing the sampling consists of correcting for the bias ^(xj 1 ) exp [— ■] . However, the variance associated to the 
unbiasing estimator in Eq. (|58[) would be too large to consider estimates reliable^ 8 .. This well-known fact results from 
the lack of overlap between the sampled and measured distributions. We therefore carry out a series of simulations 
for a set of a values ranging from to a maximum value a max so to ensure overlap between successive sampled 
distributions. 

Our choice is then to use the multistate Bennett acceptance ratio (MBAR) instead of the estimator of Eq. (|58|). 
MBAR is a method elaborated by Shirt and Choder a 48 i 58 , which aims at minimizing the statistical variance associated 
to the estimates. Following these authors, we briefly expose the principles of this procedure in a context of biased 
path ensembles. 

To use the MBAR method in the waste-recycling framework of Sec. IIVB1 we take as probability weights corre- 
sponding to a given a ensemble at each Markov chain step m the marginal probability PmargiC™) 01 Eq. (|5o]) . This 
is possible because, thanks to the detailed balance of Eq. (|17]). P^arg is preserved. 

Once one knows weights exp [— S a ] related (via Eqs. ([56)) and (|42j) ) to Pmarg((a) f° r each a-ensemble, reactivity 
averages (A) , for every ensemble a' ^ a can be computed resorting to the importance sampling identity 

(A exp [S a >]) a = Z^ ( 
(Aexp[-S a }) a , Z a [0J) 

where we used the partition functions Z a = J 2?£exp [— >S a (C)]- For a set of K different values of the bias a, a set 
(namely, a Markov chain) of M. a buffer trajectories are sampled for each bias value. An estimate of C a (t) is given by 
the MBAR estimator K for the waste-recycling averaged reactivity A a of Eq. ([55jl as 

M a 

[-4] = E W m ,cJZ (60) 

in— 1 

where weights W m ,a are given by the expression 

w m , a = z-' ^ . . (01) 



exp [-£%] 



and Z a are estimators for the partition functions Z a with minimal asymptotic covariance (see Eq. (??) and Ref. 58 ). 
Note that the denominator in Eq. (|61[) indicates that each weight W m<a takes into account contributions from all 
other ensembles k — 0, a, K. 

Introducing also partition functions Zj^ a = f T>£Ae~xjp [—S a (£)], the uncertainty can be estimated as 

var(K Q L4]) w K Q [yl] 2 (var(Z^J + var(Z Q ) - 2cov(Z^ Q , Z Q )). (62) 

Canonical equilibrium values of the correlation function Co(t) and corresponding values of the reaction rate constants 
can be recovered once we consider estimates in Eq. (|62p with a = 0. We emphasize that this is by now the first 
application of the MBAR unbiasing method based on marginal probabilities, able to estimate observables computed 
with a path sampling algorithm supplemented by waste-recycling. Numerical recipes to obtain these estimations have 
been furnished by J. Choderaff. 
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V. NUMERICAL RESULTS 
A. LJ 38 

In order to test the Lyapunov biased TPS algorithm (LyTPS), wc first consider a Lennard- Jones 38 cluster. This 
well-known benchmark system is often uysed to assess the efficiency of sampling algorithms. It and has been widely 
explored in literature for its rich thermodynamic properties: its potential energy landscape presents two main basins: a 
deep and narrow funnel containing the global energy minimum, a face-centered cubic truncated octahedron structure 
(FCC), and a separate, wider, funnel leading to a large number of icosahedral structures (ICO) of slightly higher 
energies. Although the configuration with the lowest potential energy corresponds to the FCC one, the greater 
configurational entropy associated with a large number of local minima in the icosahedral funnel make this second 
configuration much more stable at higher temperatures. As temperature increases, LJ38 undergoes several structural 
transitions. First, a solid-solid transition occurs at T ss — 0.12^- when the octahedral FCC structure gives place to 
the icosahedral ones. Secondly, above T s i = 0.18-^, the outer layer of the cluster melts, while the core remains of 
icosahedral structured 

The Lennard- Jones potential is given by the usual expression 



i (KK , v) "X 

3<k 



. 12 

er \ a 



r jk J \ r jk 



(63) 



where q J = {qi,qi,q 3 z ) is the position of the j-th atom, rjk = | q- 5 — q fe | is the distance between atoms j and k, e 
is the pair well depth and 2 1//6 cr is the equilibrium pair separation. In addition, a trapping potential prevents the 
evaporation of clusters particles at finite temperature: if the distance between the position q of a particle and the 
center of the trap q c exceeds 2.25c, this particle feels a potential |q — q c | 3 - 
The bond-orientational order parameters Qi defined a a 30 i 31 



1/2 

^ Ylm {Qjki&jk) 



Ql \ 21 + 1 N? ^ 

m— — l 



(64) 



where the Yi m (9, (f>) are spherical harmonics and 8jk,4>jk are the polar and azimuthal angles of a vectors corresponding 
to bonds between the Nf, pairs of atoms. The Qi are traditionally used to identify different crystalline symmetries. 
The parameter Q4 is often used to distinguish between the icosahedral and cubic structures, for which it has values 
around 0.02 and 0.18 respectively d 

Monte Carlo sampling fails to equilibrate the two funnels, and global optimization methods are unable to find its 
global energy minimum 1 . Hence, several elaborated algorithm have been employed in the past to study the ther- 
modynamic equilibrium of this system, such as parallel tempering^ - — basin-sampling techniques, 20 Wang-Landau 
approaches^ or path-sampling methods! 11 ^ 22 ' 23 

Standard transition path sampling 22 and discrete path sampling 24 (DPS) have been already used to study transitions 
between the two funnels of LJ38. However, in the case of TPS the large number of metastable states separating the 
two main basins prevented the traditional shooting and shifting algorithm to identify reactive paths, despite previous 
success for smaller LJ clusters^ Authors had to resort to a two-ended approach linking the two minima to find 
trajectories with the same energy of those found by DPS approach*^ The main drawbacks of this TPS method were 
a lack of ergodicity and a very large computational cost. 

Conversely, DPS has been more succesfull in this task. This method uses eigenvector following and graph transfor- 
mation— to compute the overall transition rate between two regions of phase space. To the best of our knowledge, this 
is by now the most successful approach to computing reaction rates in hJ^g,'^ In particular, reaction rate constants 
for transitions between the two solid structures have been computed using DPS 24-26 at different temperatures. 

We use here the Lyapunov biased TPS algorithm to investigate structural transitions in LJ38 for temperatures 
above and below the solid-solid transition temperature T ss — 0.12, spanning a temperature range from T = 0.10 to 
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T = 0.15. Our simulations required about 10 2 hours of cpu time to observe reactive trajectories between the two main 
funnels. Reaction constants, computed with the method exposed in Sec. lIVl can be compared to values obtained with 
the discrete path sampling approach. 22 

1. FCC-ICO reactive paths 

In order to thermalize the system at a given temperature, Langevin dynamics are run for 1000 time steps, using a 
friction parameter j5t = 1. The last configuration is used as the starting point of the first deterministic trajectory of 
the TPS simulations. At each new temperature, a new preliminary Langevin dynamics is performed. 

In simulations, trajectories consist of Af = 700 time steps, each step of duration St = 10~ 2 . Deterministic trajectories 
are obtained with the Verlet algorithm 5 ^, and then selected following the Lyapunov biased TPS algorithm described in 
Sec. IIIII For each temperature, 25 different biased path distributions are sampled, for values of the control parameter 
a ranging linearly from a = 1 to a = 2500, in order to obtain reactive paths for a = 2500 and ensure a sufficient 
overlap between distributions sampled for different a values. The unbiased distribution corresponding to a = has 
been simulated with TPS as well. A Markov sequence of 5000 biased TPS shooting and shifting moves is performed, 
in order to ensure an ergodic sampling. 

Values for the control parameters are chosen after observing the magnitude of the Lyapunov indicators L(z) for few 
trajectories, and the difference between Lyapunov indicators L(z) for current and trial trajectories in the shooting 
step (Eq. (1551 ). in order to have an acceptance ratio for the shooting move not below 20%, see Fig. [TJ The choice of 
the trajectory length r depends not on the Lyapunov bias, but on the necessity of having long enough trajectories to 
link the two funnels, r needs also to be long enough to recover an appropriate statistics to compute reaction contants 
from correlation functions (see below). The use of the Stoltz algorithm (Eq. ([30])) in the shooting moves ensures 
that the energy distribution imposed by the preliminary MD is maintained along the simulation. The value for e in 
Eq. (|30|) is taken as 0.95, so to have decorrelated sampled paths and ensure a sufficient acceptance ratio, see Fig. [2] 

We focus on the octahedral to icosahedral (FCC-ICO) transition at temperatures from T = 0.10 to T = 0.15: 
observing this passage using a direct MD or a standard TPS would require a considerable amount of CPU time 
(about 10 5 h, see Ref.— ) as the FCC configuration is at low temperatures the most stable one, so that the system 
rarely escapes from the FCC basin. In contrast, with our biased TPS technique we were able to observe at T = 0.12 
the first FCC-ICO reactive trajectories after about 300 Markov chain steps. 

To ensure that reactive paths start in the stable FCC state, we include in the path probability weight the constraining 
function (see Eq. ([25"])1 

^ cc (x ) = exp {-| (Q 4 (x ) - Q! cc ) 2 } (65) 

assigned to the starting state Xo of the path, function of the bond order parameter Q4 and centered on the value 
qFCG _ 01% We set k — 500, a sufficiently small stiffness that lets the trajectory starting point span the whole 
FCC basin. The function (p FCC keeps the beginning of the trajectories inside the FCC funnel, thus counterbalancing 
the effect of the local Lyapunov bias, that would pull trajectories on barriers. 

We present in Fig. [3] histograms for the first and the last point of the trajectories, for different values of the control 
parameter a at the FCC-ICO cohexistence temperature T = 0.12. As a values increase, trajectories explore regions 
that are increasingly distant from the initial FCC basin, and some of them eventually cross the transition region and 
reach the ICO basin. 

Once reaction paths have been identified, the computation of the inter-funnel reaction constant by the correlation 
function via Eq. (|50j) of Sec. lIVl is possible if reactants and products basins are adjacent, i.e. if there is no intermediate 
state between therm 45 ' 57 However, this hypothesis is not valid for the FCC-ICO transition: several results reported 
in the literatur o 14 ' 17 show that reactive paths linking FCC and ICO states pass through many short-lived metastable 
basins, separated by barriers of different heights, not belonging to the two main funnels. These metastable states and 
transition regions have also been observed in a previous work using the transition current sampling method/^ 2 - Such 
a feature has been confirmed as well by an attentive analysis of our trajectories. 
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Among all the intermediate metastable states, we emphasize the presence of a basin related to a faulted FCC 
configuration, indicated with D in the following, having a bond order parameter value around Q 4 = 0.12, already 
acknowledged in precedent studies ! 12 ' 17 i 72 This basin has a rather important occupation probability if compared to 
other metastable states, and is visited by every reaction path linking FCC to ICO state. Moreover, this metastable 
state is indeed structurally related to the FCC basin, and the barrier separating the D structure from FCC is lower 
than the one separating the former from ICO state. As a result, several recrossing events of trajectories starting in 
FCC, visiting the D state and then going back to FCC, can be observed. 

Hence, in order to correctly reconstruct the FCC to ICO transition paths, we have to take into account this 
intermediate metabasin. We therefore split the FCC-ICO passage in two steps: the first part is given by the passage 
from the FCC basin to the D basin corresponding to Q4 — 0.12. The second part is then given by trajectories starting 
from the D basin, and ending up in the ICO funnel. 

To obtain this second part of FCC-ICO reactive paths, we constrain the first point of the trajectories to start in 
the D basin, using a constraining function given by an indicator on the bond-order parameter value Q^xq): 

, , fl 0.10 <Q 4 < 0.13 , 
h d (Q 4 ) = { (66) 
[ elsewhere 

In Fig. 2] we present histograms for the distribution of the beginning and the end point of trajectories constrained 
with the indicator function of Eq. (|66|). at a temperature T = 0.13 slightly above the solid-solid cohexistence. Paths 
sampled with the unbiased distribution a = completely remain in the "window" given by hd{Qi(^o))- On the 
contrary, trajectories weighted with a Lyapunov bias tend to leave the metabasin: their starting points Xo tend to 
accumulate on the borders of the region defined by the indicator function in Eq. (|66|) . while the end points x T fall 
both in the FCC and the ICO funnels. 

In simulations performed at lower temperatures, this reconstruction of the second part of the FCC-ICO reactive 
path with trajectories starting from the D state is more difficult. Indeed, histograms for trajectories of the same 
length at temperatures lower than the solid-solid transition T = 0.12 show that an important fraction of the sampled 
trajectories fall from the D state directly to the FCC basin, while a few trajectories end in the ICO state. This is 
attributed to the heights of the barriers separating the metastable D structure from either the stable ICO or stable 
FCC structures, the latter barrier being lower than the former one. 



2. FCC-ICO reaction constants 



The total reaction constants for the two-step FCC-ICO transition, assuming a steady occupation probability for 
the intermediate D state, is derived in Appendix [X] and reads 

k F ^ d kd^i 

kp ^ = 71 TT \ ( 67 ) 

where subscripts F, d and I refers to FCC, D and ICO states respectively. The same steady state approximation is 
assumed for all intermediate metastable states in discrete path sampling studies* 1 - ' 24 ' 28 

Reaction rates kF^-d and kd^i involve transitions between states separated by high free energy barriers,— thus the 
hypothesis of time scale separation required by the reaction rate theory is still valid, and reaction constants can be 
computed using the method exposed in Sec. IIV1 Reactive paths between FCC and D basins, and between D and ICO 
basins, are computed as reported above (Sec. IV AT]) . On the contrary, the D to FCC reaction rate kd-yF cannot be 
computed by LyTPS, because the requirement of a time scale separation is no longer valid, the barrier separating this 
two states being too low. It is therefore computed by direct MD simulation. 

The reactivity A (Eq. (fSU)) ) for each trajectory is evaluated in simulations distinguishing the three basins FCC, ICO 
and D whose ranges of bond-order parameter Q 4 value, that is 0.13 < Q± < 0.18, < Q4 < 0.04 and 0.1 < Q4 < 0.13, 
respectively. Data harvested during LyTPS runs are unbiased using MBAR. 
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T 


kF^d 






k F ^i (Ref. 1 ) 


0.10 


1.2 10~ 7 


1.4 10~ 7 


8.1 10~ 14 


2.5 10~ 13 


0.11 


1.3 10~ 5 


2.5 10~ 7 


1.08 10" 11 


1.15 lO -11 


0.12 


8.1 10~ 5 


4.0 10~ 7 


1.2 10~ 10 


2.82 10~ 10 


0.13 


3.3 10~ 4 


2.0 10~ 5 


6.6 10~ 9 


4.2 10 -9 


0.14 


9.3 10" 4 


4.5 10" 5 


4.3 10~ 8 


4.3 10 -8 


0.15 


2.4 10 -3 


2.4 10 -4 


5.7 10~ 7 


3.2 10~ 7 



TABLE I: Table of reaction constants for the transitions from FCC to D structure, D to ICO, and the total FCC to ICO 
transition, indicated as kp^d, fcd->/ and kF^i respectively, at different temperatures. Values of fc_F->/ are obtained using 
Eq. (|A2|) and assuming the reaction constants fcd-s-F as 10 -2 , 3 • 10~ 2 , 10~ 4 for T — 0.10, T = 0.11 and T — 0.12 respectively 
(values obtained by Langevin MD), and unitary for T > 0.12. In the last column on the right, we report DPS data from Ref.—, 
computed in the harmonic approximation framework. 

In Fig. [5] and [6j two examples of population correlation functions for the computation of reaction rate constants, 
unbiased with MBAR, are reported. Note that reactivity values computed at short times are nearly zero, and do not 
contribute significatively to the correlation functions: in fact, these values are obtained from segments of trajectories 
too short to witness a complete transition between two states. In Table QJ we report reaction rate constants values 
for the FCC to D structure (kF^d) transition, and fro the D structure to ICO (kd-+i), that give, through Eq. (|A2I) . a 
total FCC to ICO (kp->i) rate in good agreement with values given by DPS calculations 24,28 . Finally, an Arrhenius 
plot comparing our results with the reaction constants proposed in Refill is presented in Fig. [7] 

B. Vacancy migration in a-Iron crystal 

The second example of a thermally activated process studied using LyTPS is the migration of a single vacancy 
in a-Iron crystal. Atomic interactions of the model system are described by an embedded atom model potential^! 
The crystal structure is body-centered cubic, and the initial unrelaxed cell contains 1023 atoms displayed on 1024 
lattice sites, the vacant site corresponding to the vacancy. The reaction coordinate used to represent the motion of 
the vacancy is the distance crossed by the moving atom that replaces the vacancy. 

The free energy landscape for this system has been reconstructed in Ref.—. It presents two stable states, the first 
corresponding to the initial configuration, and the second one to the same configuration obtained by translating a 
neighboring atom into the vacancy site. The first neighbor distance is a = 2.47A, switching its initial position with 
the vacancy site. There is a single free energy barrier separating these two states for temperatures above T — 450K, 
while for lower temperatures an intermediate metastable state appears, corresponding to an intra-site position for the 
moving atom«2& 

We performed LyTPS simulations with trajectories of different lengths (see below), with time step 8t = 4- 10~ 15 s. 
A preliminary MD simulation is done to equilibrate the system to the required temperature, with a friction parameter 
7 = 2.5 • 10 12 s~ 1 . We explored temperatures ranging from 300^ to 850i\\ The LyTPS shooting and shifting moves 
are iterated to cosntruct Markov chains consisting of 1500 trajectories. 

As for LJ38, the trajectory length and the values for the control parameter a have been chosen in order to ensure 
an acceptance ratio of 25% and an adequate ergodic sampling of the phase space. For temperatures above 450i^, 
the presence of a single "smooth" barrier separating the two metastable states makes this application simple enough: 
sampling of reactive trajectories is achieved using 15 a values from the unbiased simulation at a — up to a — 150T0 12 , 
with trajectories of 300 steps. For temperatures below 450-fT, an ergodic sampling of trajectory space appears more 
difficult. We therefore employed longer trajectories of 500 time steps, as well as larger values of the control parameter, 
up to a = 500 • 10 12 to allow the system to escape the initial basin. 

Reaction constants for the passage between the two stable states above 450-ftf are estimated from correlation 
functions unbiased with the MBAR algorithm, via Eq. (|50]). 

For T < A50K, the presence of an intermediate metastable basin has to be taken into account in the evaluation 
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of reaction constants. As recalled in Sec. IV Al the reaction rate expression obtained from Eq. (|50|) holds only for 
adjacent reactant and product basins. At low temperatures, it is therefore more appropriate to use our algorithm to 
evaluate the reaction constant for the passage from the initial state to the intermediate basin. From Eq. (|A2[) , reaction 
rate for the passage from one stable configuration to the other is simply half the reaction constant from one stable 
configuration to the intermediate one. Indeed, reaction constants for transitions from the intermediate metastablc 
state to either of the two stable states are equal, because of the symmetric shape of the potential surface.— 

In Fig. [8] we compare the reaction rates obtained with LyTPS, those computed inserting in the transition state 
theory (TST) expression (Eq. ([52]) ) the free energy barriers reported in ReL 7 - 3 -, and reaction constants obtained with 
a classical harmonic approximation (HA). Above the Debye temperature (A70K), rates obtained with LyTPS fall 
between TST and harmonic approximation values. To explain this point, we recall that reaction rates we estimate 
with the method exposed in Sec. IIVI are derived from Eq. (|50[) . hence correspond to the phenomenological rate 
constants. These values are therefore bounded from above by TST values, that overestimate reaction rates^i, as can 
be seen from Eq. (|51[) . Conversely, values obtained with the harmonic approximation neglect anharmonicity effects on 
the activation barrier, thus giving reaction rates that are lower then the phenomenological rate constants we compute. 
Our results are then in agreement with the reaction rate theory recalled in Sec. IIVI 

VI. CONCLUSION 

The method presented in this paper allows to compute reaction rate constants for inter funnel transitions in many- 
body systems. The reaction rate values are evaluated using a path sampling algorithm biased with local Lyapunov 
numbers. This bias aims at enhancing an accelerate sampling of reactive paths, so as to reduce the lenght of Markov 
chains and the amount of CPU time to observe activated processes. We assess the performace of the method by 
observing reaction paths and evaluating equilibrium reaction rates for structural transitions in the LJ38 system and 
for vacancy migration in an a-Iron crystal. For both systems, we were able to predict phenomenological rate constants, 
in very good agreement with data already given in the literature in the case of LJ38. 

The Lyapunov biased TPS method presents several advantages, and incorporates features of different rare events 
simulation methods. 

Firstly, with respect to other importance sampling methods based on Lyapunov weighted sampling 8 ' 33 , Lyapunov 
biased TPS has the main advantage of a simpler implementation. This is due to the Lyapunov indicator L(z) we 
propose in Eq. (|22l) . that allows to quantify the divergence of hamiltonian trajectories by resorting to local Lyapunov 
numbers. These quantities can be easily calculated with the Lanczos algorithm, that enables one to compute the 
largest eigenvalues of the Jacobian matrix of the hamiltonian mapping with a limited computational cost. As recalled 
in Sec. A|) . resorting to local Lyapunov numbers to evaluate chaoticity of phase space trajectories doesn't suffer from 
the computational drawbacks of other algorithms aimed at the same purpose, as RLI or the tangent space method.— 
The implementation of shooting and shifting Monte Carlo moves in a Lyapunov biased TPS is less complicated, and 
computationally less expensive, than the algorithm proposed in Ref. 33 with the use of RLI, because we do not need 
to compute four trajectories to evaluate the divergence of a single pat h 33 ' 43 . 

Secondly, this formulation for the Lyapunov indicator is such that the bias applied to each path in order to 
enhance the fraction of reactive trajectories is clearly identified, contrarily to what is done in rather complex cloning 
algorithms like the one proposed in Lyapunov weighted dynamics^ and transition current sampling 7 ^. Hence, the use 
of standard unbiasing statistical tools to recover unbiased observables is possible and requires a small theoretical and 
computational effort. 

Furthermore, we consider the access to the evaluation of equilibrium transition rates as the most important aspect 
of Lyapunov biased TPS. On the computational point of view, the direct access to reaction rates without resorting to 
a distinct evaluation of the reaction barriers and the transmission factor, as usually done in standard TPS technique 2 , 
is a very advantageous feature. To unbias reaction constants computed in Lyapunov biased ensembles we chose 
among other unbiasing algorithm, like WHAMZ£ or Extended Bridge Sampling 76 techniques, the MBAR method 48 . 
MBAR has proven to be computationally efficient and to give an adequate numerical precision in estimating reaction 
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constants. Moreover, this work is the first in which MBAR is implemented exploiting the marginal probability derived 
from a waste recycling method. This allows to include informations encoded in rejected proposal trajectories. 

Finally, LyTPS is implemented at a finite temperature, imposed to trajectories by the canonical distribution from 
which the path starting point is selected and maintained along the path thanks to the Stoltz proposal for the shooting 
algorithm, see Sec. lIIII In parallel, the Lyapunov indicator used as a bias to select reactive paths directly links the path 
sampling to the local conformation of the potential energy surface via the hessian matrix, thus giving to our method 
an intrinsic dependence on the potential energy landscape. The coupling between a finite-temperature sampling and 
potential energy surface conformation is a noticeable improvement if compared to eigenvector-following methods, that 
are based on the shape of the potential energy surface, but usually operate at zero temperature. LyTPS can be 
acknowledged as a finite, nonzero temperature version of the well-known eigenvector-following techniques, such as 
Dimer, Optim or ART^l^ 

The advantages related to a finite temperature technique do not concern only the exploration of the energy land- 
scape, but also the fact that the evaluation of physical observables like reaction rates takes into account temperature 
and anharmonicity effects. Indeed, the phenomenological reaction rate we computed (see Sec. IIV|) can be compared 
with experimental measures: LyTPS turns out to be a powerful tool to study problems like vacancy migration in Iron. 
In this context, reaction rates are usually estimated using only the potential energy barriers at K and harmonic 
approximations: this give approximate results with respect to experimental data obtained at nonzero temperatures. 
LyTPS reaction rates find an important application as input parameters for object Monte Carlo codes aimed at 
numerically reproducing resistivity recovery experiments. 

We conclude observing that this method can be furhter developed using a parallel tempering technique, or imple- 
mented for the computation of reaction rates in more complex condensed matter systems, and can find interesting 
applications in a wide class of research fields, spanning from molecular biophysics to physical metallurgy, where the 
numerical determination of reaction rates has important consequences for experimental applications. 
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Appendix A: Reaction constants for processes with intermediate states 



The time evolution of the occupation probabilities p F (t), Pi(t) an d Pd{t) of the FCC, ICO and faulted states 



respectively, read 



p F (t) = -k F ^ d p F (t) + k d ^ F p d (t) 

Pd(t) = k F ^ d p F (t) - (fc rf ^F + kd^-i)Pd(t) 

pi(t) = k d ^!p d (t) - fc/^ d p/(£) 



ki^ d pi{t) 



(Al) 



where subscripts F, d and / refers to FCC, faulted and ICO states respectively, and fc/n-s indicates the generic 
transition rate from state A to state B. Note that the possibility of a direct transition from FCC to ICO states 
without passing by the intermediate basin has been neglected. The stationary approximation for the metastable state 
imposes p d (t) = 0, thus the above system can be recasted in the form 




- tf > Pf(*) + t? M lV F \ Pi(t) 

(kd^F + k d ^l) rr v 1 (kd^F+kd^i)^ 1 v > 
(kd^F+kd^ri)^ r y 1 (k d ^F+kd^l)^ ly ' 



and the reaction constant between FCC and ICO states can be calculated as 

k F ^ d k d ^i 



k F . 



(k 



d^F 



(A2) 
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FIG. 1: Top: Average value of the Lyapunov indicator L(z) over a Markov chain of 1000 trajectories, starting from the FCC 
basin at T = 0.12, as a function of the control parameter a used in the simulations, for different trajectory lengths r. Increasing 
a increases the mean Lyapunov indicator and enables trajectories to explore barriers and transition states. Average Lyapunov 
indicators are almost independent of the trajectory length r. Bottom: Acceptance ratio, given by Eq. (|35p for the same 
simulations. 
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FIG. 3: Top: Histogram of the initial point position xo for trajectories starting from the FCC basin at T = 0.12 for different 
values of the control parameter a, averaged on a Markov chain of 5000 steps. The restraining function of Eq. (|65[) mantains 
the initial states of the trajectories in the FCC funnel for all a values. Bottom: Same histogram, for the final position x T . 
Trajectories sampled with large a values escape the FCC funnel more often. Their final states are distributed over the whole 
FCC-ICO range. 
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FIG. 4: Top: Histogram of the initial configuration xo for trajectories starting from the D configuration metabasin located 
at T — 0.13 for different values of the control parameter a, averaged over a Markov chain of 5000 steps. Bottom: Same 
histogram, for the final position x T . Trajectories end up in both the FCC or the ICO funnel. 
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FIG. 5: Left: Correlation function for the transition from FCC to D basin, at T = 0.13. Error bars are computed via Eq. (|62[) 
and directly given by Chodera's numerical recipe^. Right: Reaction constant for this same passage, obtained at times shorter 
than the mean first passage time. The reactive flux k(t) reaches a plateau value, corresponding to fcF->d, as explained in Sec. lIVI 
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FIG. 6: Top: Correlation function for the D to ICO transition, at T = 0.13 . Bottom: Reaction constant for this same 
passage. 
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FIG. 7: Arrhenius plot for the FCC to ICO reaction rate from Table |T](LyTPS, red dots) compared with data obtained from 
DPS 24 using an harmonic approximation and reported in Ref.— (HA, black line). 
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FIG. 8: Arrhenius plot of reaction constants for migration of monovacancy in a-Iron obtained with Lyapunov biased TPS 
(LyTPS, red points), compared with rates obtained using in Eq. (|52[) the free energy barriers proposed in Ref.— (/ctst, black 
line) and using an harmonic approximation (HA, green line). 



